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Abstract 

In this paper, we present a new hybrid algorithm for the time integration of col- 
hsional A^-body systems. In this algorithm, gravitational force between two particles 
is divided into short-range and long-range terms, using a distance-dependent cutoff 
function. The long-range interaction is calculated using the tree algorithm and inte- 
grated with the constant-timestep leapfrog integrator. The short-range term is calcu- 
lated directly and integrated with the high-order Hermite scheme. We can reduce the 
calculation cost per orbital period from 0{N'^) to O(A^logA^), without significantly 
increasing the long-term integration error. The results of our test simulations show 
that close encounters are integrated accurately. Long-term errors of the total energy 
shows random-walk behaviour, because it is dominated by the error caused by tree 
approximation. 

Key words: methods: n-body simulations — solar system: formation 
1. Introduction 

For the time integration of collisional A^-body systems such as star clusters and systems 
of planetesimals, the combination of direct summation for force calculation and the individ- 
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ual timestep algorithm has been the standard method for nearly half century (Aarseth 1963, 
2003). For galactic dynamics and cosmology, fast and approximate methods for force calculation 
such as the particle-mesh scheme (Hockney & Eastwood 1981), the P^M scheme (Hockney & 
Eastwood 1981), the tree method (Barnes & Hut 1986) and combinations of PM and tree (Xu 
1995, Bagla 2002, Dubinski et al. 2004, Springel 2005, Yoshikawa & Fukushige 2005, Ishiyama 
et al. 2009) are used. It is not impossible to combine individual timestep algorithm and fast and 
approximate force calculation. For example, McMillan & Aarseth (1993) developed a high-order 
integrator using individual timestep combined with the tree algorithm. However, it was difficult 
to achieve good performance on distributed-memory parallel computers for such scheme. Fujii 
et al. (2007) introduced a hybrid of tree and individual timestep algorithm, which is designed to 
handle the evolution of star clusters embedded in the parent galaxy. In their BRIDGE scheme, 
both the parent galaxy and the star cluster are expressed as A^-body systems. The interaction 
between particles in the star cluster are calculated directly and integrated with the individual 
timestep scheme, while interactions between particles in the parent galaxy and that between 
particles in parent galaxy and particles in star clusters are calculated with the tree algorithm 
and integrated with the leapfrog scheme with shared and constant timesteps. 

The BRIDGE scheme is based on the idea of splitting the Hamiltonian of an A^-body 
system to multiple components. The time integration of tree part is symplectic and does not 
generate any secular error. The direct integration of the internal motion of star cluster is not 
symplectic, but treated with high accuracy using high-order integrators. The obvious limitation 
of the BRIDGE scheme is that it can handle the close encounters of particles in star clusters 
only. If we want to allow close encounter between particles in the parent galaxy, it goes back 
to the usual direct summation scheme. 

For the time integration of planetary systems, the mixed variable symplectic (hereafter 
MVS) scheme (Kinoshita, Yoshida & Nakai 1991, Wisdom & Holman 1991) has become the 
standard method. For the time integration of almost stable orbits of planets, the MVS is 
well suited. However, if we want to handle protoplanets or planetesimals, we need to handle 
their close encounters and collisions. The MVS scheme, however, cannot handle them since 
it requires that the timestep is kept constant (Calvo & Sanz-Serna 1993). In order to handle 
close encounters, several modifications of MVS method have been proposed. One is the hybrid 
method by Chambers (1999). It handles close encounters with high accuracy by splitting the 
Hamiltonian to that of close interaction and distant interaction as in the case of P^M method. 
It can be used for planetary accretion problems. However, it relies on the direct calculation 
and its calculation cost scales as 0(A^^). Very recently, Moore, Quillen & Edgar (2008) applied 
GPGPU to the hybrid method. In addition, they used the Hermite scheme (Makino & Aarseth 
1992) for the integration of the short range force. 

Brunini & Viturro (2003) and Brunini et al. (2007) combined the hybrid method and 
the tree algorithm. Therefore their method can handle a large number of particles and close 
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encounters. The scheme described in Brunini & Viturro (2003) uses the leapfrog integrator 
to handle the Hamiltonian for the short-range interaction. Therefore it had to use very small 
timesteps. The scheme described in Brunini et al. (2007) seems to be improved, but no details 
of the scheme was given. 

In this paper, we describe a new time integration algorithm which combines the strong 
points of these schemes. This scheme is combination of the BRIDGE scheme (Fujii et al. 2007) 
and a hybrid symplectic integrator method (Chambers 1999). We call this scheme as Particle- 
Particle Particle- Tree (hereafter PPPT). In section 2 we overview previous numerical methods. 
These methods are based on the MVS method. Therefore, we first introduce the MVS method 
and then overview extension of MVS method. In section 3 we describe our new time integration 
algorithm, PPPT. In section 4, we show the result of test simulations. In this paper, we discuss 
only the simulations of planet formation process. However, in future works, we apply PPPT to 
other collisional systems. Summary and discussions are given in section 5. 

2. The Numerical Method 

2.1. The symplectic integrator 

The Hamiltonian of an A^-body system is given by 

N 2 N^l N r-mm- 

j 2mj .^-^ j=i+i 

where pt is the momentum of particle i, rij is the distance between particles i and j, mi is the 
mass of particle i and G is the gravitational constant. The Hamilton's equation of motion is 

(2) 

where {f,H} is the Poisson bracket and / is a canonical variable. We define the differential 
operator as Df = {f,H}. The general solution of equation (2) at time t + At is formally written 

as 

/(t + At)=e^*^/(t). (3) 

Equation (3) cannot be solved analytically. We can obtain the approximate solution by dividing 
the Hamiltonian into multiple parts that can be analytically solved. 

In the symplectic integrator, the Hamiltonian is divided into two parts as 

H = Ha + Hb, (4) 

i=i j=i+i 

(6) 

where Ha is the potential energy and Hb is the kinetic energy, and each of which can be solved. 
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The time evolution of / is given by 

/(t + At)=e^*^/(t)=e^*(^+^)/(t), (7) 

where A = {,Ha} and B = {,Hb} are operators. The exponential in equation (7) can be 
approximated as 

k 

^AtiA+B) ^ -Qga,AMgfe.AtB ^ 0{Ar+^), (8) 
i=l 

where (ai, 02, a/c) and (61, 62, are real numbers, k is the number of stages, n is the 

order of approximation. The first term of the right hand side of equation (8) is an order-n 
approximation of the left hand side. The first-order symplectic integrator is given by 

/(t + At)=e^*^e^*^/(t)+0(At2), (9) 
and the second-order symplectic integrator, which is the leapfrog integrator, is given by 

f{t + At) = e^*^/2e^*^e^*^/VW + 0(At3). (10) 

The symplectic integrator has the advantage that there is no long-term energy error for the time 
integration of periodic systems. In the case of a near-Kepler potential, both the semi-major axis 
and the eccentricity are conserved. One disadvantage of the symplectic method is that high- 
order schemes are expensive and have rather large local error coefficients. Another disadvantage 
is that it requires a constant timestep. If we change the timestep following usual recipes for 
Runge-Kutta-Fehlberg-type schemes (Fehlberg 1968), the energy is no longer conserved (Skeel 
& Gear 1992, Calvo & Sanz-Serna 1993). There are a number of proposals for methods to 
combine the variable timesteps and the good nature of the symplectic schemes. Most of them 
are based on the idea of splitting the potential to fast-varying and second-varying terms and 
applying different timesteps. For example, Skeel & Biesiadecki (1994) proposed a method that 
splits the gravitational force into many components, each of which has finite effective range in 
the distance. By assigning different timesteps to different components, they effectively realized 
a variable timestep integration. 

The methods discussed below are all based on the idea of changing the way to split the 
Hamiltonian. 

2.2. The Mixed Variable Symplectic Method 

Kinoshita, Yoshida & Nakai (1991) and then Wisdom & Holman (1991) introduced the 
MVS method for planetary systems. In the symplectic scheme, we split the Hamiltonian to 
kinetic and potential energy, so that we have analytical solutions for both parts. This division 
is not unique. As far as each splitted Hamiltonian has an analytic solution, any division can be 
used. The idea of MVS is to split the Hamiltonian to the Keplarian term H^ep and interaction 
term Hjnt- The time evolution is given by 

H = Hxep + HjnU (11) 



fit + At) = e^*^/2g^*7gAii^/2^^^)^ (^2) 

where K is the operator defined as Kf = {f,HKep} and / is the operator defined as // = 
{f,Hjnt}. The MVS method integrates Hjnt with the leapfrog integrator and Hxep by using 
analytic solution of the Kepler orbit. Thus the MVS method is expressed as follows: 

1. Calculate accelerations due to gravitational interactions between planets at time t and 
give a velocity kick. 

2. Update analytically positions and velocities from t to t + Athj using the Solar gravity. 

3. Calculate accelerations due to gravitational interactions between planets at time t + At 
and give a velocity kick. 

The advantage of MVS is that only interactions between planets are integrated numerically. 
The Solar gravity is analytically integrated and is accurate up to the round-off error. 

2.3. The BRIDGE code 

The BRIDGE code was introduced by Fujii et al. (2007). It was designed for the time 
integration of star clusters embedded in parent galaxies. This scheme is a combination of the 
direct and a tree schemes, using the idea similar to that of the MVS scheme. The internal inter- 
actions of stars in star clusters are integrated by the Hermite integrator with direct summation, 
while other parts are integrated by the leapfrog integrator and tree scheme. 

The BRIDGE scheme divides the Hamiltonian as 

H = H^ + Hp, (13) 

^GmciniGj ^^Gmcimcj .... 
^» = -2^— 1^1^— ' (14) 

i<j ' GG,ij j=l j=l ' GC,ij 

where Nq and Nsc are the number of particles in the parent galaxy and that in the star cluster, 
mG,i and pc.i are the mass and the momentum of particle i in the galaxy, mc j and pc,i are those 
of particle i in the star cluster, and rcGi'f'GC^fcc are distances between two galaxy particles, 
one galaxy and one star cluster particle, and two star cluster particles, respectively. We can 
express the time evolution from t to t + At as 

f{t + At) = e5^*"e^*^e5^*°/(t), (16) 

where a is the operator defined as a/ = {f,Ha} and /3 is the operator defined as /?/ = {f.Hp}. 

The BRIDGE code uses the leapfrog scheme for if^, and the fourth-order Hermite scheme 
for Hp. Thus the integration procedure during a tree timestep At is done in the following way: 

1. Make a tree at time t and calculate accelerations from all particles on galaxy particles, 
and from galaxy particles on star-cluster particles. 

2. Give a velocity kick for star cluster particles, and update the velocities of galaxy particles. 
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3. Integrate positions and velocities of star cluster particles from t to t + At by using the 
Hermite scheme with the individual timestep, and positions of galaxy particles by making 
them drift with the constant velocities. 

4. Make a new tree at t + At and calculate accelerations from all particles to galaxy particles, 
and from galaxy particles to star cluster particles. 

5. Give a velocity kick for star cluster particles and update velocities for galaxy particles. 

In this scheme, is integrated with the symplectic leapfrog scheme, while Hp is integrated 
with non-symplectic Hermite scheme. It combines the fast tree code for the orbital motion of 
particles in the galaxy and high-accuracy Hermite scheme for the internal orbital motion of 
particles in the star cluster without any additional approximation. Thus, it is the first scheme 
with which we can follow the orbital and internal evolution of a star cluster embedded in a 
galaxy in a fully self-consistent way. 

24. The MERCURY code 

The MERCURY code (Chambers 1999) splits the Hamiltonian into three parts as 

H = Hxep + Hint + Hsun, (17) 

H...-i(^-^)-i^W(r,), (18) 

i=l \^"H 'i I i^j 'ij 

Hint = -j:^^{i-W{n,)), (19) 

i=l ^"^0 

where Hjnt is the potential energy of gravitational interactions between particles except for 
those undergoing close encounters, H^ep is the kinetic energy of particles plus potential energy 
of particles undergoing close encounters, Hsun is the kinetic energy of the Sun, and W{rij) is 
the changeover function. The modification from MVS is that the potential energy of nearby 
particles is separated from Hjnt and moved to H^ep- The time evolution is described as 

fit + At) = e5^*^e^^*^e^*^^e5^*^e^^*^ /(t). (21) 

where / is the operator defined as // = {f,Hj}, S is the operator defined as Sf = {f,Hs} and 
K is the operator defined as Kf = {fjHx}- This scheme is sometimes called the hybrid scheme. 
In the actual code (MERCURY), the force from gravitational interactions is split instead of the 
Hamiltonian for the ease of programing. We can write the pairwise force F[rij) as 

F{r,j) = F,iose{rij) + F^isMj), (22) 
F,i,se{rij) = F{r,^K{rij), (23) 
FaisMj) = F{nj)[l-K{r,^)], (24) 
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where Fciose{fij) is the force for close encounter and Fdist{fij) is the remaining force. The 
relation between W and K is given by 



In the hybrid scheme, Hxep dose not have an analytical solution if the potential of close en- 
counter is not zero. Therefore, Kf is integrated numerically using the Bulirsch-Stoer method 
(Bulirsch & Stoer 1964 and Stoer & Bulirsch 1980). The one-step integration of the hybrid 
scheme is done in the following way: 

1. Apply the velocity kick, e^^*^, due to distant interaction Hint- 

2. Apply the position drift e^^*"^, due to Hsun- 

3. Integrate orbits of particles which are under close encounters using BS method, and update 
positions and velocities of the rest of particles using the Kepler orbit. 

4. Give the position drift by Hsun with stepsize At/2. 

5. Give the velocity kick by Hjnt with stepsize At/2. 

The hybrid method is used for long-term integrations of outer solar system, restricted 
three-body problems, and planetary embryos. It has a good performance for accuracy and speed 
for small- systems, but for system with a large number of particles it becomes expensive. 

2.5. The DAEDALUS code 

Brunini et al. (2007) presented a new mixed-variables symplectic tree code for plan- 
etesimal dynamics, DAEDALUS. This code is an improved version of the modified tree code 
described in Brunini & Viturro (2003). Brunini & Viturro (2003) developed the tree code 
(Barnes & Hut 1986) with two-level timesteps for planetesimal systems. This tree code usually 
integrates all particles with a constant timestep using the leapfrog integrator. If close encounters 
occur, it integrates only particles which undergo close encounters with much smaller timesteps. 
When integrating close encounters, timesteps are determined using equation (2) in Brunini & 
Viturro (2003), and the tree is constructed at each steps. 

The DAEDALUS integrator is combination of the tree code and a hybrid symplectic 
integrator method. It splits the Hamiltonian following the description of Chambers (1999). 
Therefore the DAEDALUS integrator integrates Hjnt using tree method and integrates Hxep 
analytically where there are no close encounters. If close encounters occur, it integrates Hxep 
numerically using the Bulirsch-Stoer method. The one-step integration of the DAEDALUS 
integrator is done in the following way: 

1. Make a tree and calculate accelerations from Hjnt then give a velocity kick. 

2. Give a position drift by Hsun with stepsize At/2. 

3. Calculate positions and velocities for particles which are in close encounters, and update 
positions and velocities with the Kepler orbit for particles which are not in close encounters. 

4. Apply the position drift e^^*'^'^ due to Hsun- 




(25) 
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Table 1. Characteristics of each method. 



5. Make a new tree at next step and calculate accelerations from Hjnt then give a velocity 
kick. 

The DAEDALUS integrator uses the variable, but shared, timestep for close encounters. 
Therefore, when the number of particles in close encounters is large, the calculation cost can 
become high. If we use an individual timestep, we can decrease the calculation cost significantly. 
In Table 1, we summarize the characteristics of the previous and present algorithms. Table 1 
shows that there were no algorithms which has all of the listed desirable properties except 
for the one described in this paper. In this paper, we present a new algorithm which uses 
the shared timestep for distant interactions and the individual timestep for close interactions. 
Furthermore, the force due to distant interactions are calculated by using tree method with a 
changeover function. 

3. Particle-Particle Particle- Tree (PPPT) scheme 

In this section we describe our new scheme for collisional A^-body systems. We use the 
fourth-order Hermite method for near-neighbour forces and the tree method for distant forces, 
and we use the hybrid method to split the gravitational force by using a changeover function. 
In our scheme, we split the gravitational force between two particles as 

F{U,) = FHardinj) + Fsoftinj), (26) 

FHard{rij) = F{rij)K{rij), (27) 
FsofMj) = F{r,j)[l-K{r,^% (28) 

where F^ardi'^ij) is the force from nearby particles and Fsoft{fij) is the force from distant 
particles. These formulae are the same as equations (22)-(24). We can write the Hamiltonian 

as 

H = Hnard + Hsoft, (29) 

/ i<j 




Hsoft = L — - Win,)), (31) 



where Hnard contains the kinetic energy of all particles and the potential energy of near- 
neighbours, and Hsoft contains the potential energy of all other pairs of particles. We treat the 
solar gravity as the force caused by a fixed potential in this paper. So we do not have Hsun 
here. We can express the time evolution from t to t + At as 

fit + At) = e5^*^e^*^e^^*^/(t). (32) 

The changeover function K splits the gravitational force between particles to contribu- 
tions of close encounters and others. Thus, Hnard changes rapidly and Hsoft changes slowly. In 
this paper, we use two types of changeover functions. One is the fourth order spline function 
which introduced by Abe et al. (1986) given as 

A>.)=(2^^)\ (33) 

where X = nrij/rcut and Vcut is a scaling radius. Note that K{rij) becomes zero where Vij > Tcut 
(see fig 1). The other was first introduced by Levison & Duncan (2000). It is given by 

[1 if 1" > 1, 

Kivij) = I lOY^ - 15Y^ + 6Fi° if < F < 1, (34) 

[o ify<o, 

where Y = -^^z^- Hereafter we call it the DLL function. In figure 1, r\jrcut = 0.4 and r2/rcut = 
0.6. We call r2 the cutoff radius of the DLL function. 

In this paper we regard the gravitational field of the Sun as an external potential. This 
treatment is okay for the study of planet formation process of earth-type planets, because 
planetesimals do not perturb the Sun strongly. The total mass of planetesimals is much smaller 
than the mass of the Sun and planetesimals are distributed almost uniformly around the Sun. 

We integrate Hsoft using the leapfrog and the tree method with a constant timestep At. 
We integrate HHard using the fourth-order Hermite method with the block timestep. For the 
timestep criterion, we used a slightly modified version of the "standard" criterion (Makino & 
Aarseth 1992 and Aarseth 2003). The standard criterion is given by 





II (2) 


1 + 


\di\ 


-1 


\di\\ 


(3)| 


+ 1' 


(2) 


2 



V \ai\\al '\ + \al ^ 

where rj is the accuracy parameter. The new timestep criterion we used is given as 





(2) 


1 + 1 


di\'^ 




Ml 


1 + 1 


(2) 


2 



^^^ = VJ , . ,, ,3^, . , (o),. , (36) 



ao = a — —. (37) 

Here Oq is a constant introduced to prevent the timestep from becoming unnecessarily small 
when I Oil is small, and a is a parameter which controls the size of the timestep. Here, is the 

(2) (3) 

acceleration due to Fnard, and di,a- and a- are its first, second and third time derivatives, 
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respectively, and rrii is the mass of particle. With equation (35), the timestep becomes unnec- 
essarily small if there is just one particle inside the radius Vcut of one particle and Vij ^ rcut- 
To illustrate this problem, consider the case in which one particle moves away radially with a 
constant velocity v. High order derivatives of the force from this particle is given by 

a = FK, (38) 

a = {F'K + FK')v, (39) 

^(2) _ ^^//^ ^ 2F'K' + FK"y, (40) 

^(3) _ ^^w^ ^ ^pf'j^' ^ 3^/^// ^ FK"'y (41) 

for the equation (35), where F is the gravitational force from a particle and K is the changeover 
function and F'{x) = dF/dx. In the case of < Vcut, "we can expand K around Vcut- Without 
the loss of generality, we can assume that r^ut = 1 and dr/dt = v = 1. Then we have 

K = -Z^ + 0{Z^), (42) 

K' = -5Z^ + 0{Z^), (43) 

K" = -20Z^ + O{Z^), (44) 

K'" = -60Z^ + O{Z'^), (45) 

where Z = . By substituting equations (43)-(45) into equations (39)-(41), and omitting 

time derivatives of gravitational force, we obtain 

a = -Z^F + 0{Z^), (46) 

a = -5Z'^F + 0{Z^), (47) 

a(2) = -20Z3f + O(Z^), (48) 

a^^') = -Q0Z^F + O{Z^). (49) 

Because time derivatives of gravitational force are individual to Z. Equation (35) becomes 



AbZ^F^ + OjZ^) 

^ 700Z6F2 + O(Z7)' ^ ^ 



rjZ ^ 



9 + 0(Z) 



\\ 140 + O(Z)' 

= r]Z{12Q0 + O{Z)). (52) 

Equation (52) shows that the timestep approaches to zero as r^j approaches to Vcut- This 
behaviour is clearly undesirable, since there is no need to reduce the timestep for the neighbour 
force, when the neighbour force itself is small. The reason why criterion (35) gives zero stepsize 
is, as we can see from equations (47)- (49), \a\ approaches to zero faster than its high order 
derivatives. However, this is due to the cutoff by the changeover function, and the actual 
physical acceleration by the particle just inside radius Vcut is of the order ^J^. In order to avoid 
this unnecessarily small timestep, we introduce in equation (37). By doing so, we make 
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criterion (37) to give the timestep which is accurate relative to the absolute strength of the 
force itself, before the changeover function is applied. This timestep criterion does not use the 
gravitational force from the Sun. In other words, the timestep is determined purely by the 
forces from nearby particles. If we included the Solar gravity, the original criterion could lead 
to unnecessarily small timesteps, since the Solar gravity is much larger than the forces from 
neighbour particles. 

Our scheme is summarized as follows: 

1. Make a tree at time t and calculate accelerations due to Hsoft- 

2. Give a velocity kick. 

3. Integrate positions and velocities from t tot + At using the Hermite scheme with the block 
timestep and Hfjard- 

4. Go back to step 1. 

When we integrate Hnard, we use the list of neighbours for particles to save the calculation 
time. If we do not use the neighbour lists, we have to calculate forces from all particles and 
the calculation cost becomes 0{N'^). We construct the neighbour list of particle i by selecting 
particles within distance r„/ from particle i at time t. Here Tni must be sufficiently larger than 
Tcut , SO that particles outside the radius r.„; do not enter the sphere of radius Vcut during one 
timestep. In this paper, we use r„; ^ Vcut + 3Atcr, where a is the the velocity dispersion. In 
order to find neighbours fast, we use a uniform 2-D grid with the grid size smaller than Vcut- To 
summarize our neighbour finding way, we first assign all particles to cells. We then look over 
the neighbouring cells which are within r„; from a particle. The particles in the neighbouring 
cells are its neighbours. 

4. Accuracy and Performance 

In this section, we present the result of test calculations for the accuracy and performance 
of our new algorithm. We adopted the distribution of planetesimals following the Hayashi model 
for test calculations. The surface density at 1 AU is 10 g/cm^. The unit mass, the gravitational 
constant G, and the unit length are normalized to one solar mass, one, and 1 AU respectively. 
For most of the tests, we use 10000 equal-mass particles distributed randomly between the radii 
of 0.9 and 1.0 AU. Their mass is 1.45 x 10^^ g and their velocities follow the Rayleigh distribution 
with <e>= 5rH, <i>= 2.5r/^, where < e > and <i> are the dispersions of the eccentricity and 
inclination. Initial radius of particles is 364 km, which corresponds to the density of 3 g/cm^. 
Physical collisions are handled under the assumption of perfect accretion. The radius of the 
collision product is determined to keep the density unchanged. Unless specified otherwise, all 
test calculations are for 10 orbital periods. 

Figures 2 and 3 show the relative energy error of the system, \E — Eq\/\Eq\ where 
is the energy of the system at time 0, and the number of direct gravitational interactions per 
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one tree timestep per one particle, Ndirect, as a function of the cutoff radius Vcut, for the case 
of the sphne changeover function (eq. 33). The cutoff radius Vcut is normahzed by the Hill 
radius th at 1 AU. In figure 2, we plot the largest energy error during the time integration 
for lOyr (10 periods). We use t] = 0.05, timestep At = 0.040 — 0.0050 yr, and opening angle 
6 = 0.5 and 0.1. The relative energy error is practically independent of Vcut if ^cut/^H > 3 and 
At < 0.020 yr. The number of mutual interactions is proportional to the square of the cutoff 
radius (figure 3). Since the scale height of the disk is about brn, it is in most cases smaller than 
rcut- Therefore the number of particles is proportional to the square of radius. Figures 2 and 3 
show that by using Vcut/^H ~ 3, we can achieve high accuracy with very small value of Ndirect- 
On the other hand, the increase in the calculation cost is pretty small even when we use very 
large Tcut/^H such as 50, because the force calculation using the tree is more expensive. In the 
case of 6' = 0.5, the energy error is lower-bounded at 10~*, while with 9 = 0.1, the error can be 
reduced to lO^^''. In this case, using At = 0.040 yr resulted in a rather large error as shown 
in figure 2. This behaviour can be understood by looking at the time evolution of the error. 
Figure 4 shows the time evolution of the error for At = 0.040 yr and 0.010 yr. The behaviour 
of the error shows quasi-periodic behaviour with a period ~ 1 yr for At = 0.040 yr. In figure 2, 
the relative energy error for At = 0.040 yr, Tcut/f^H = 20 is larger than that for At = 0.040 yr, 
fcut/'^H = 7. This is because the random energy error from tree scheme with large timestep is 
dominant for At = 0.040 yr. Thus, this peculiar behaviour occurred. 

Figures 5 and 6 show the energy error for the case of the DLL cutoff functions, with the 
inner radius ri/rn = 1 and 10. We can see that the behaviour of the error is quite similar to 
that in the case of the spline cutoff, and independent of the choice of ri. 

Figures 7 and 8 show the relative energy error and the number of interactions per one 
particle in the tree part as a function of the opening angle 6. In figure 7, the energy error 
shows the power-law dependence as oc 6'^'^ for the case of At = 0.0050 yr. On the other hand, 
for At = 0.040 yr, the error dose not go below 10~^ for small values of 6. This is because, 
for At = 0.040 yr, the truncation error of the integrator becomes larger than the error due 
to force approximation for ^ = 0.1. In figure 8, the number of interactions per one particle is 
about two orders of magnitudes smaller than that for direct calculation, for ^ = 0.1. In figure 
8, we can see that the dependence of Nint to 6 is rather weak. In the case of stellar systems, 
calculation cost is proportional to 6*"^ ^ (Makino 1991). In our experiment the dependence is 
~ 6*"^. This is because the distribution of particles is a thin and narrow ring. Figure 9 shows 
the relative energy error as a function of opening angle for the case of the DLL function again, 
the behaviour of the error is similar to the case of spline function. 

Figure 10 shows the relative energy error as a function of the size of timestep for Hsoft- In 
the case of 6* = 0.1 and At < 0.010 yr, the error is dominated by that from the tree approximation, 
and becomes independent of At. If we want to keep the error in lOyr (10 periods) to be less 
than 10~^, a pair of 6' ~ 0.2 and At ~ 0.02 is probably a good choice. In realistic calculations 
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10~ in 10 orbits is probably okay, though how small the error should be is a difficult question. 
At least, the error of order of 10~^ is smaller than that of energy change of planetesimals 
due to gas drag and coUisional damping. We do not show the result for the DLL cutoff here. 
We calculated by using same parameters for the DLL cutoff and conffimed that the result is 
essentially the same as that for the spline cutoff in figure 10. 

Figures 11 and 12 show the relative energy error and the calculation cost of neighbour 
force Ndirect as functions of the accuracy parameter rj. Other parameters are 6 = 0.1, rcut/rn = 10 
and At = 0.0050 yr. Figure 11 shows that the energy error is practically independent of the 
choice of a. On the other hand, in figure 12, small a results in the increase of Ndirect- In 
practice, a pair of r/ = 0.1 and a = 1 seems to be a good choice. 

Figure 13 shows the long-term variation of the relative energy error. The calculation is 
done with the opening angle 9 = 0.5, the cutoff radius Tcut/^H = 10 and the timestep At = 0.0050 
yr. In this case, the energy error reaches about 7.5 x 10~^ after the time integration for 10^ yr 
(10^ orbital periods), while it is 9.8 x 10~^ for 10 yr (10 periods). In other words, the energy 
error grows 10 times larger as the integration time becomes 1000 times longer. It shows that 
the growth of energy error is stochastic like a random walk. It means that the error is mainly 
caused by the force error of the tree scheme (Barnes & Hut 1989). The growth of energy error 
is, therefore, expected to be slow and the error is small enough even after long calculations. 

Figure 14 shows the calculation time per one tree timestep as a function of the total num- 
ber of particles in the system N. The calculation is done with the opening angle = 1, the cutoff 
radius rcut/fn = 5 and the timestep At = 0.0050 yr. We used the Intel(R) Core(TM)2 Quad 
CPU Q6600(2.4GHz). It shows that the calculation time increases as O(A^logA^). Therefore, 
we reduce the calculation cost from 0{N'^) to O(A^logA^). 

Figure 15 shows the number of tree interactions, Nfree per particle as a function of A^. 
We can see that Ntree is roughly proportional to O(logA^). 

5. Summary and Discussion 

We have developed a new hybrid A^-body simulation algorithm for the simulation of 
coUisional A^-body systems. This new scheme is constructed by combining the tree and direct 
schemes using the hybrid integrator. The results of test simulations of evolution of a planetesi- 
mal system show that our new scheme PPPT can drastically reduce the calculation cost, to the 
level comparable to the cost of a tree scheme with constant timestep while keeping accuracy 
sufficient for realistic simulations. 

In principle, our scheme can be used for coUisional systems other than planetary systems, 
such as globular clusters or stars around a supermassive blackhole in the galactic center. We'll 
show the results of simulations of such systems using our scheme in future. 

We are grateful to Tomoaki Ishiyama, Masaki Iwasawa, Michiko Fujii, Kuniaki Koike, 
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Keigo Nitadori and Yusuke Tsukamoto for fruitful discussions. This work was partially sup- 
ported by the Research Fund for Students (2009) of the Department of Astronomical Science, 
the Graduate University for Advanced Studies, KAKENHI 21244020 and 21220001. 
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Fig. 1. The changeover functions. The horizontal axis is the distance between i and j particles normalized 
by the cutoff radius, and K(rij) is the changeover function. The solid and dashed curve, show the fourth 
order spline function and the DLL function. This function uses ri/rcut = 0.4 and r2lrcut = 0.6. 
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Fig. 2. The energy error plotted against the cutoff radius. Crosses, triangles, squares and circles show 
the results with At = 0.04,0.02,0.01 and 0.005 yr, respectively. The left and right panels show the results 
with = 0.1 and 0.5, respectively. 
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Fig. 3. The number of direct gravitational interactions per one tree timestep per one particle with the 
spline function plotted against the cutoff radius. Crosses, triangles, squares and circles show the results 
with e = 0.1, A< = 0.04,0.02,0.01 and 0.005 yr, respectively. 
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Fig. 4. The relative energy error of the system with the spline function plotted against the calculation 
time. The solid curve shows tlic result with — 0.04 yr, t^.^i /ru~7 and = 0.1. The dashed curve shows 
the result with At = 0.04 yr, rcutj^H = 20 and Q = 0.1. The dotted curve shows the result with At = 0.01 
yr, rcut/rn = 7 and 9 = 0.1. 
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Fig. 5. The relative energy error of the system with the DLL function plotted against the r2 cutoff radius. 
Crosses, triangles, squares and circles show the results with At = 0.04,0.02,0.01 and 0.005 yr, respectively. 
The left and right panels show the results with ti/th = l,?'2/''ff = 2 — 100,0 = 0.1 and 10,11 — 100,0.1, 
respectively. 




Fig. 6. The relative energy error of the system with the DLL function plotted against the r2 cutoff radius. 
Crosses, triangles, squares and circles show the results with At = 0.04,0.02,0.01 and 0.005 yr. respectively. 
The left and right panels show the results with ri/rH — l,r2/rH — 2 — 100,0 = 0.5 and 10,11 — 100,0.5, 
respectively. 
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Fig. 7. The relative energy error of the system with the sphne function plotted against the opening angle 
6. Open squares, open circles, filled squares and filled circles show the results with rcut/rn = 10 and 
At = 0.04 yr, rcut/rn = 50 and At = 0.04 yr, rcut/rn = 10 and At = 0.005 yr and rcut/rn = 50 and 
At — 0.005 yr, respectively. 




e 

Fig. 8. The number of tree gravitational interactions per one tree timcstep per one particle with the 
spline function plotted against the opening angle 0. Open squares, open circles, filled squares and filled 
circles show the results with r cut/ I'm = 10 and At 0.04 yr, rcut/rn = 50 and At = 0.04 yr, Vcut/rn = 10 
and At = 0.005 yr and rcut/rn = 50 and At = 0.005 yr, respectively. The four results are practically 
indistinguishable. 
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Fig. 9. The relative energy error of the system with the DLL function plotted against the opening angle 
9. The left panel shows the results with ri/ru = 1 and open squares, open circles, filled squares and filled 
circles show the results with r^jrH = 3 and A< = 0.04 yr, r2/rH = 10 and At = 0.04 yr, ^2/^/^ = 3 and 
At = 0.005 yr and r2/r^f = 10 and At = 0.005 yr, respectively. The right panel shows the results with 
Tih'H = 10 and open squares, open circles, filled squares and filled circles show the results with r2/rH = 15 
and At = 0.04 yr, r2/rH = 50 and At = 0.04 yr, r2/rH = 15 and At = 0.005 yr and r2/rH = 50 and 
At = 0.005 yr, respectively. 




At [yr] 

Fig. 10. The relative energy error of the system with the spline function plotted against the tree timestep. 
Open squares, open circles, filled squares and filled circles show the results with 9 = 0.1 and rcut/rn = 10, 
9 = 0.1 and rcut/rn = 50, 9 = 0.5 and rcut/rn ~10, 9 = 0.5 and rcut/i"H = 50, respectively. 
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Fig. 11. The relative energy error of tlie system plotted against the timestep accuracy parameter rj. 
Crosses, squares and triangles show the results with a = 1,0.1 and 0.01, respectively. 
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Fig. 12. The number of direct gravitational interactions per one tree timestep per one particle with the 
spline function plotted against the timestep accuracy parameter rj. Crosses, squares and triangles show 
the results with a = 1,0.1 and 0.01, respectively. 
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Fig. 13. Energy error of the system plotted against a function of time. 
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Fig. 15. The number of tree gravitational interactions per one tree timestep per one particle with the 
spline function plotted against a function of number of particles. Crosses, open squares and open triangles 
show the results with 9 = 0.1,0.5 and 1, respectively. 
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